## Heteroskedastic Probit via ML in R

  hetprob<-function(X,Z,y, method='BFGS'){
    X<-cbind(1,X)

    neglnl<-function(theta,X,Z,y){
        b<-theta[1:ncol(X)]
        g<-theta[ncol(X)+1:ncol(Z)+ncol(X)+1]
         print(b)
         print(g)  
	p<-as.vector(pnorm((X%*%theta)/exp(Z%*%g)^2))
	-sum(y*log(p) + (1-y)*log(1-p))
        }        
        
    result<-optim(c(mean(y),rep(0,ncol(X)-1),rep(0,ncol(Z))),
      neglnl, hessian=T, method=method, X=X, Z=Z, y=y)

    coef<-result$par
    var<-solve(result$hessian)
    se<-sqrt(diag(var))
    zscore<-coef/se
    pvalue<-2.0*(1.0-pnorm(abs(zscore)))

    print(cbind(coef,se,zscore,pvalue))
#    print(var)
    print(result$convergence==0)
    }
  

